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1.0 SUMMARY 
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An analysis technique for simulation of supersonic mixed- 
compression inlets with large flow field perturbations (hammershock, 
un start/ restart, etc.) is presented. The approach is based upon a 
quasi-one-dimensional inviscid unsteady formulation which includes 
engineering models of unstart/ restart, bleed, bypass, and geometry 
effects. Numerical solution of the governing time-dependent equations 
of motion is accomplished through a shock-capturing finite-difference 
algorithm, of which five separate approaches are evaluated. Comparison 
with experimental supersonic wind tunnel data is presented to verify 
the present approach for a wide range of transient inlet flow 
conditions. 
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2.0 INTRODUCTION 


The function of a supersonic inlet is to supply the airflow required 
by an engine at the highest possible pressure level while maintaining 
minimum drag. In order to minimize iniet cowl drag for sustained flight 
at speeds above Mach 2.0, it becomes essential that some portion of the 
supersonic area contraction be accomplished internally. An inlet of this 
type is commonly referred to as a mixed-compression inlet. For mixed- 
compression inlets, optimum internal performance is provided by main- 
taining the terminal shock at the inlet throat. This operation provides 
high-pressure recovery and minimizes distortion at the engine face. 
These inlets, however, have a discontinuous airflow characteristic 
known as unstart. If an airflow transient causes the terminal shock to 
move upstream from the throat, the shock is unstable and is abruptly 
expelled ahead of the cowling. This shock expulsion, or unstart, 
causes a sharp reduction in mass flow and pressure recovery and a 
large drag increase. Inlet buzz, compressor stall, and/or combustor 
blowout may also occur. Obviously, an inlet unstart is extremely 
undesirable because of the adverse effects not only on the propulsion 
system itself but also on the aerodynamic qualities of the aircraft. If 
inlet unstart does occur, complex mechanical variations to alter the inlet 
geometry are required to restart the inlet and stabilize the normal 
shock at a position downstream of the throat. 

Prior modeling efforts of supersonic inlet gas-dynamic phenomena 
have concentrated on a control volume or lumped-parameter approach 
following the. original work of Martin (Ref. 1). This type of simulation 
is basically a one-dimensional mathematical model (involving ordinary 
differential equations) for predicting the transient behavior (dynamic 
response) of the subsonic duct downstream of the terminal normal 
shock, with the subsonic duct divided into a number of lumped 
volumes. Models that treat the whole subsonic duct as a single lumped 
volume may be sufficiently accurate for fairly short air induction 
systems, but their accuracy deteriorates for inlets which are long. For 
such long air induction systems, it becomes necessary to divide the 
subsonic duct into more than one lumped volume. Although simulation 
accuracy increases with increasing number of volume lumps used, 
practical considerations limit this number to three or four. 




The supersonic inlet investigation by Amin and Hall (Ref. 2) 
divided the subsonic duct volume into three initially equal volumes. 
The volume nearest to the throat had a moving upstream boundary (the 
normal shock) and a fixed downstream boundary; the remaining two 
lumped volumes had fixed boundaries. The instantaneous total tempera- 
ture and total pressure within each of these volumes were obtained 
during a simulation run by numerically integrating the corresponding 
rate equations. These rate equations were derived by applying the 
one-dimensional unsteady continuity and energy equations in conjunction 
with the equation of state to each lumped volume. Complete derivation 
of these equations is given in Appendix 1 of Ref. 2. To calculate the 
rate of change of mass flow at each of the lumped volume boundaries, 
the rate of change of momentum within each volume due to the instan- 
taneous net imbalanced force acting or' it during transient conditions 
was utilized. Complete derivation of the mass flow rate equations is 
given in Appendix II of Ref. 2. For the case of inlet unstart, the 
continuity method by Moeckel (Ref. 3) was used to calculate the 
expelled shock position ahead of the inlet cowl and the relationship 
between this shock position and the resultant subsonic flow spillage 
over the inlet cowl. 

A linearized characteristics type of analysis has been developed by 
Willoh (Ref. 4) for dynamic response studies of one-dimensional inviscid 
inlet flows. This approach combines a set of linearized equations across 
the normal shock with an exact solution of the linearized wave equation. 
In general, the Willoh analysis is more exact than conventional 
lumped-parameter techniques and, for many problems, is no more compli- 
cated in application. However, it is based upon a linearized treatment 
and hence is strictly applicable only to inlet flows with small flow-field 
perturbations. Results presented in Refs. 5 and 6 show reasonable 
agreement with experiment for mixed-compression inlet frequency 
response (amplitude ratio and phase angle) at a free-stream Mach 
number of approximately 2.5. 

The work by Mays (Ref. 7) was one of the first to employ a 
completely numerical solution of the one-dimensional, unsteady, inviscid 
flow equations in a variable area duct. His numerical technique was 
based upon the Lax artificial viscosity explicit algorithm. Through this 
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approach, large amplitude transients (such as compressor surge) and 
their tuect on mixed-compression inlet flow could be numerically simu- 
lated. The main difficulty with this approach is the use of the Lax 
algorithm which is well-known to introduce spurious second-order 
dissipative (viscous-like) terms into the governing equations. However, 
for the time in which this work was performed (1968), the Lax algorithm 
was the best numerical technique available. 

The present report will document an analytical investigation into 
simulation of large perturbation flow field effects (such as hammer- 
shock, un start/ restart, etc.) in mixed-compression supersonic inlets. 
A quasi-one-dimensional inviscid unsteady approach will be described 
which includes engineering models of unstart/ restart, bleed, bypass, 
and geometry effects. Numerical solution of the governing time- 
dependent equations of motion will be accomplished through a shock- 
capturing finite-difference algorithm, of which five separate approaches 
are evaluated. Comparison with experimental supersonic wind tunnel 
data will be presented to verify the present approach for a wide range 
of transient inlet flow conditions. 

The following advisory panel provided high-level technical exper- 
tise and guidance during the course of the present work. 
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Fluid Dynamics 
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Inlet Engineering 
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3.0 ANALYSIS 


An illustration of an axisymmetric mixed-compression supersonic 
inlet is presented in Fig. 1. Pertinent details which will be modeled 
by the present work include: 

• cowl bleed and bypass 
centerbody bleed 

• centerbody translation 

• supersonic inviscid flow field between the tip and the inlet 

plane 

• variable boundary conditions at the exit plane 

The X- and y-axis will be taken as shown (along and normal to the 
centerline of the centerbody, respectively) . 

3.1 GOVERNING EQUATIONS 

If the cross-sectional area of a flow passage varies very slowly 
and the radius of curvature of the central axis of the passage is large 
contrasted to the passage height, the flow inside the passage is said to 
be a quasi-one-dimensional flow. In this case, the flow area ^ is a 
function of both distance x and time t. The flow properties are 
assumed to be uniform across all surfaces perpendicular to the mean 
flow direction. Figure 2 illustrates the quasi-one-dimensional flow model 
for unsteady inviscid flow with friction fJ, heat transfer ' q , and mass 
bleed mj. 

The derivation of the appropriate governing fluid dynamic equa- 
tions of motion (continuity, momentum, and energy) may be found in 
most compressible fluid flow textbooks, for example. Chapter 24 of 
Shapiro (Ref. 8) and Chapter 19 of Zucrow and Hoffmann (Ref. 9). 

These equations are (with reference to Fig. 2 for nomenclature): 
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FIGURE 2. FLOW MODEL FOR UNSTEADY OUASI-ONE-DlMENSIONAL INVISCID FLOW 
WITH FRICTION. HEAT TRANSFER. AND MASS ADDITION 


Continuity 

iiM), + a M 
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Momentum 

at" ax ^ ax ^ 

Energy 

+ iLMliiil = 3A ^ ^ 
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' iiT*) * p(c^T ♦ iiT*) 

= Mass bleed term 

= Friction term ♦ mass bleed momentum term 

= Heat transfer term ♦ mass bleed energy term 
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The flow is taken to be thermally and calorically perfect air obeying the 
perfect gas equation of state 

P = - (y - 1) [e - iTfd*] (4) 


with 


Y 
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For the present supersonic inlet flow application, friction and heat 
transfer effects are not considered. 

Equations (1), (2), and (3) are written in what is known as weak 
conservation law form, with source/sink terms on the right-hand side. 
When these equations are written in this form and are solved numeric- 
ally for a supersonic flow, embedded shock waves and expansion waves 
can form and decay automatically without special treatment of any kind. 

The source/sink terms M , F , and Q are used to incorporate mass 

s s s 

bleed/bypass effects as well as inlet un start/ restart phenomena. 
Further observe that the area variation terms in the momentum and 
energy equations also play the role of source/sink-like terms. 


3.2 NONDIMENSIONAL GOVERNING EQUATIONS 

The governing equations of motion (1), (2), a.id (3) may be 
written in the following vector form 

iT+ir=‘r (6) 

3t 3X 

where the vector components U, F, and G are given by 
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G = pA^/A + Fg 


rpA^/A + Q, 
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with the following nondimensional variables 




X 


1 s *^s*~ref 
^ ^ref^ref^ref 


^ = — £| — 
^ref\ef 


A u 


^ref^refVef 


t = EA 

‘’ref^ref^ef 


'ref 


ft = *^s*‘ref 
^ ‘^ref^ref^ref 


F = *^s*~ref 
s •^ref'^ref^ref 


In the above, terms denoted by "ref" denote constant reference 
conditions, selected to be free-stream sonic conditions in the present 
work. In addition, the reference area is related to the reference 
length through 

A s L * 

^ref ^ref 

and the equation of state (4) is given by 

p = (y- 1 ) [e - ip5'] no) 

in terms of the nondimensional variables. 

The remainder of this report will utilize subscript notation to 
denote partial differentiation. In subscript notation, Eq. (6) becomes 


lit ♦ F, ■ 6 


. 1 . 












where 
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3.3 COORDINATE TRANSFORMATION 

The governing vector equation (11) is now transformed in terms of 
the generalized coordinates t and C where 


$ = 5(x, t) (12) 

For a function f = f(x, t) 



where the subscript denotes partial differentiation, i.e., 

f . 8f f - ii 

“ 35 * t“ 3t 


Since = 0 from the definition of x, the Jacobian matrix of the 
generalized coordinate transformation becomes 



(14) 


In terms of the generalized x, C coordinates, the governing vector 
equation (6) may be written as 


(x^u)^ * (5, x^ F ♦ . x;r 

Recalling that x^ = J from Eq. (14) and defining 

U = x^= J U 
F = 

6 = 

the transformed governing vector equation (15) becomes 
“t* ''x*" * ® 

12 


(15) 


(16) 













(17) 
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where the term is known as the metric coefficient and the term K as 
the grid speed (see Ref. 10 for the clarification). ^ 

The above transformation of coordinates can be thought of as a 
real nonsingular mapping from the physical (x, t) space to a computa- 
tional (5, t) space. Equation (17) will be solved by a numerical method 
(to be described later) in the computational (?, x) space with the 
resulting solution transformed back into physical (x, t) space for 
interpretation. For the present work, the grid speed is taken to be 
zero; thus, the resulting 5 grid distribution will not be allowed to 
change with respect to time, i.e., the grid is nonadaptive. A good 
discussion of adaptive grids and the many problems involved in proper 
determination of the grid speed may be found in Ref. 10. 

Under the restriction of zero grid speed, Eq. (17) reduces to 

^ (18) 

where 

X J (19) 

is related to the Jacobian of the transformation given in Eq. (14). For 
implementation in the various numerical algorithms to be discussed in 
the next section, a grid spacing of unity (AC=1) is utilized in the com- 
putational space. The Jacobian J in Eq. (19) is then computed numeri- 
cally per the following schematic illustration. 

M - “ 1 . - 1 , 


M J 1+1 

COMPUTATIONAL SPACE 




AX - X, - X,_, 

AX = Xj+,-*, 

lu] 

1 


1-1 


j+1 


PHYSICAL SPACE 
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3.4 NUMERICAL METHODS 


3.4.1 Beam-Warming Algorithm 

The Beam-Warming solution technique is based on the work of 
Beam and Warming (Refs. 11 and 12). The governing flow equations as 
given in Eq. (18) are solved using a non-iterative, implicit, first- 
order-time-accurate solution technique constructed in delta format (i.e., 
increments of the conserved variables and fluxes). The delta format is 
used extensively to achieve analytical simplicity and numerical efficiency 
in addition to the advantageous property of a steady state, when one 
exists, independent of the time step. Fourth order dissipation is added 
explicitly to the solution variables to achieve numerical stability over a 
wide range of Courant numbers. 

The following subsections describe the Beam-Warming technique as 
applied to the solution of the one-dimensional flow equations as given in 
Eq. (18). The subsections are devoted to the definition of: 

1. Implicit solution technique 

2. Boundary point technique 

3. Boundary condition constraints 

4. Hybrid solution technique 

3. 4. 1.1 Implicit Solution Technique 

Following Ref. 12, the Fade’ formula is applied to the vector form 
of Eq. (18). The Fade' formula for implicit time differencing is 


3t 


,N 


1_ UN 

At V 1+A / 


( 20 ) 


where the Fade’ formula has been restricted to the Euler implicit 
temporal differencing with truncation errors 0(t*). This specific form 
was chosen over others available from this formula (such as the Euler 
explicit, trapazodial, and three point backward) based on a combination 
of ease of implementation, relative speed, and accuracy. As used here, 
the superscript N denotes the time level, and A is the forward time 
difference operator, given by 


AU = . uN 
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Applying the Fade' formula given in Eq. (20) to the governing 
one-dimensional flow equations, Eq. (18), there results 

+ At (F^ - = 0 (22) 

where F = C F . 

* 


The terms F and G are locally linearized in time by using truncated 
Taylor series expansions 


pN+1 . ?N „N+1 . c * X aN .. 

r “ ^ + F At + • • * = A U 


N+1 


+ F_^At + 

T 


6*^+1 = + gJ AU^ + G_At + . . . » + B 


aU^ + G At + ••• 

T 


where A is the Jacobian matrix of the governing one-dimensional equa- 
tions Eq. (18). Here the homogenous property of the flux vector was 
used to simplify the expansion of the flux vector, F. Substitution of 
the expansions given by Eq. (23) into Eq. (22) and collecting terms, 
Eq. (22) can be approximated as 


[I + At ( 


3C 


- bN)] .u" . pN 


(24) 


where 

= At[G^At - (aV)^ + G^] 

As used here, I is the identity matrix and 
operator. 




s a matrix 


(25) 


Equations (24) and (25) can be evaluated once the spatial deriva- 
tive operators are replaced by finite-difference operators. The solution 
vector at the new" time-step is obtained by casting Eqs. (24) and (25) 
into a block tri-diagonal structure and solving for the deltas of the 
dependent variables where the old and new time-step solutions are 
related to the deltas through Eq. (21). Applying central, second order 
spatial differences to Eqs. (24) and (25), Eq. (26) results 


S iUj.j 4 B + f AUj^j . 


(26) 
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where A 5 - A, B s I - AtB and C s A . Here, j can take 
on values insido the upstream to downstream boundaries. Explicit 
fourth-order dissipation is added to the right hand side, r^, terms to 
provide numerical stability following Ref. 12. These dissipative terms 
are of higher order than the difference forms used in the governing 
equations and consequently do not disrupt the formal accuracy of the 
method. Boundary point formulations similar to those for interior point! 
will be described in the following subsection. 

3. 4. 1.2 Boundary Point Technique 

The boundary point-difference equations are formulated to have 
the same time accuracy as the interior point equations and are com- 
patible with the latter in their spatial order of accuracy. This 
treatment of boundary points is patterned after that presented in 
Ref. 13 in that it provides a fully implicit set of equations for the 
boundary points. The time differencing for boundary points is the 
same as that for interior points. Hr,vever, it is convenient to perform 
the spatial differencing before linearizing the equations in time. From 
this approach, the boundary conditions enter naturally into the dif- 
ference equations. 

Associated with the outflow boundary point is a half-cell where the 
centroid of the half-cell is indicated by the asterisk as shown below. 





HALF CELL CENTROID 


EXIT PLANE NODE 
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The exit plane node point which lies in the exit plane boundary is 
designated as point J. The difference form of the governing Eq. (18) 
at the boundary half-cell centroid is 


iU* 

3t 


+ ’fj • Sj 


C27) 


where 7 is the backward spatial difference operator. Equation (27) 
retains the global conservation property as pointed out in Ref. 13 if the 
half-cell centroid is defined by linear interpolation of data between the 
centroid of the last interior cell and the boundary as 

u* ■ (1 - 5> “j (28) 

Combining Eqs. (27), (28), and the Fade' formula, Eq. (20), and 
writing the result in a block tri-diagonal structure, Eq. (29) results. 

A + B aUj = rj (29) 

with A = 1/4 “ AtA, B - 31/4 At(A“B) and rj is an Interpolated value 
using points J and J-1 . 

A similar analysis for the inflow or entrance plane boundary equa- 
tion results in the inflow boundary point equation given in Eq. (30). 

B aUi + C a'Jz * r* (30) 

Here, B = 1/4 ♦ AtA and C = 31/4 - At(A + B) and r*. is obtained by 
interpoiation between points 1 and 2 at the upstream boundary. 

The inflow and . outflow equations given by Eqs. (29) and (30) 
complete the block tri-diagonal equation structure for the boundary as 
well as interior grid points necessary to solve for the interior flow 
solution. The following subsections will consider the required outflow 
or exit plane and inflow or inlet plane boundary constraints required to 
close the system of solution equations. 
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3. 4. 1.3 Boundary Condition Constraints 

For purely supersonic outflow, all characteristics (u S a) have 
positive slope and eminate from the interior toward the boundary. The 
flow variables at the outflow boundary are determined by the interior 
flow. In this case, the full set of flow boundary conditions, Eq. (29), 
are used to evaluate boundary conditions at the outflow boundary. 

For subsonic outflow (u S a) all except one characteristic has a 
positive slope as shown below. 


QUAUli'V 



Information propagates upstream along the characteristics from the 
boundary point to the interior flow region. For the ''•ase of u S a at 
the outflow boundary, one boundary condition or constraint may be 
specified. In this case, Eq. (29) is modified to account for the con- 
straint where a linearized algebraic boundary condition is used in place 
of one of the three equations. This algebraic relation will be referred 
to hereafter as a constraint. 

Permissible exit plane boundary condition constraints, as used 
here for u S a, are specified pressure, Mach number, and mass flow. 
Corrected mass flow can be related directly to the exit plane Mach 
number and, thus, is treated like an exit plane Mach number 
constraint. 

With the aid of Eq. (4) and the definition of the U vector in 
Eq. (7), the boundary condition with static pressure p as a constraint 
will be derived. Using Eq. (4) as the algebraic relation to impose the 
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pressure p as a constraint and writing in terms of the U vector ele- 
ments, one obtains 


P = (Y-l)(U3-iU|/Ux) 


(31) 


This constraint is to be satisfied at x 


N+1 


by a function F such that 




(32) 


where ( 33 ) 

F(U.t) = (y-1)(U3-1U|/Ux) - p(t) = 0 

If the constraint function F is linearized in time and the constraint 
defined in Eq. (32) is imposed, the linearized form of the constraint 
equation can be cast in a block tri-diagonal form that is compatible with 
the governing equation structure. The result is 


(0)4Uj.j + B iUj - rj 


(34) 


N N 

with B = -|jj- ^U,t), and Tj = -F (U,t) - . This algebraic 

equation is used to replace the third equation (energy equation) in 
Eq. (29) for the exit plane boundary condition. 


The outflow Mach number constraint is imposed in a manner similar 
to the exit pressure constraint developed in the previous paragraph. 
The outflow Mach number constraint is imposed through the Mach number 
definition as given in Eq. (35) 


M = 



(35) 


and is written in terms of the U vector elements as 


F(U,2) = (Uj/Ux)/[y(y-1)(U3/Ui - ifU^/Ux)]^ - M(t) (36) 
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The constraint equation given in Eq. (34) is imposed with, however, 
the right hand side term, rj, given by 

rj , t) - 1^ (U. t) ^ IT (37) 

This algebraic constraint equation for Mach number replaces the third 
or energy equation in Eq. (29). 

For the outflow mass flow constraint, the second equation in 
Eq. (29) (momentum equation) is replaced. The mass flow constraint is 
given in Eq. (38) . 

F^(U, t) » ^2 - m(x) “ 0 (38) 

where the rj term in Eq. (34) takes on the form 

r, . -f"(U. t) - ^ (U. t) ^ iT (39) 

The inflow boundary constraints are formulated in a manner similar 
to the outflow constraints as developed in the previous paragraphs. 
Permissible boundary conditions at the inflow boundary were developed 
by examining the characteristics at the inlet plane. With reference to 
the sketch below for supersonic inflow where u > a, all characteristics 
have positive slope and are directed downstream into the interior. 
Since no information propagates upstream across the inflow boundary, 
all inflow boundary variables must be specified and are not constrained 
by the interior flow solution. 
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3.4.2 Bearri'Warming Hybrid Algorithm 

Typical supersonic inlet flows are characterized by moving 
moderate-strength shocks. For a numerical solution scheme to be 
adequate and robust as applied to inlet flows, it must resolve these 
discontinuities accurately with minimal overshoot and remain stable as 
shock strengths increase. The Beam and Warming method developed in 
the previous subsections, even with significant fourth-order dissipation 
added, can become numerically unstable for moving moderate-strength 
shocks. This is due to the characteristic pressure undershoot in the 
predicted pressure just upstream of the shock in a supersonic inlet 
flow. Increasing the fourth order dissipation will not resolve the 
stability problem and will seriously degrade the shock tracking accuracy 
of the scheme. 

In order to partially resolve this numerical resolution and stability 
problem, the spatial difference operators as applied in the Beam- 
Warming method are repiaced by switching operators that convert from 
central to upwind (one-sided) across shock discontinuities. Basically, 
this modification attempts, in a heuristic manner, to replace the central 
difference operators in the Beam-Warming technique with operators that 
reflect the propagation of information based on the local characteristic 
directions at a point in the flow. This modification to the Beam- 
Warming technique, hereafter referred to as the Beam-Warming Hybrid 
technique, follows the development given in Ref. 11. 

The central spatial difference operator used to develop Eq. (26) is 
replaced by a switching difference operator given in Eq. (40) where the 
switching or transition operators that change from central to upwind 
maintain strict conservation and local consistency. 

ly =: i_ V + ^t,/2 - 7eA/2 (40) 

95 A5 1 - 7e/2 

As used here, V and A are the ciassical backward and forward 
finite-difference operators, respectively, and e is a switching 
parameter. When e = 0, Eq, (40) reduces to a central difference 
operator, and when b = 1, an upwind difference operator results. 


Applying the Eq. (40) difference operator to Eqs. (24) and (25) 
and combining terms, Eqs. (41) and (42) result 

[7^ + (Bj_j + 

+ [I - Gj/2 + Arcj (bJ + A^)/2 + Atej_iAj/2] aU^ 

— -1 — N ^ 

■t [C “ ^^j+1 ° 


r" s r” + Itej [-gJ + (Aj..iUj,.i-a5Uj)]/2 (42) 


Here, the matrices A, B, and C, and the right hand side vector, r^, 
are the Beam-Warming matrices as used in Eq. (26). 

In a like manner, the forth order dissipation difference term in the 

4 

Beam-Warming scheme, 6 U., is replaced by a conservative switching 
operator of the form 

E-h Qj 

A 3Uj - M3) 

A 

The switching scheme is implemented by defining s = 0 for subsonic 
points where a central difference is to be imposed and e = 1 for 
supersonic points where a upwind scheme is applied. Conservation of 
this scheme is maintained across a shock discontinuity by redefining the 
switching operator for the last supersonic point upstream of the shock 
at c * 0. 


3.4.3 Flux Vector Splitting Algorithm 

The Beam-Warming hybrid solution technique outlined in the pre- 
vious subsection provides an engineering approach to account for the 
propagation characteristics of flow information in flows containing 
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moving shock discontinuities, It will be shown in a later section that 
this approach yields robust solutions for flows containing moving 
shocks. Numerical .olution oscillations still persist, however, in 
subsonic regions downstream of the shock that cannot be adequately 
damped without seriously degrading solution accuracy. 

A new solution technique recently introduced by Steger and 
Warming (Ref. 14) and extended to a shock capturing algorithm by 
Reklis and Thomas (Ref. 15) provides a physically and mathematically 
consistent approach to include characteristic regions of influence in the 
solution algorithm. In this approach the flux vector of the governing 
one-dimensional equations is split into flux pieccr corresponding to the 
three eigenvalues of the linear system. Contributions of the flux 
pieces, depending upon the signs of the eigenvalues, are employed to 
structure the matrix coefficients. A, *B, and C (similar to those used in 
Eq. 26) so that spatial differences also reflect the direction of 
information propagation. This splitting and recombination of flux 
information is entirely consistent with the propagation of characteristics 
information. It results in a method that not only yields robust solu- 
tions for moderate strength moving shocks but also exhibits enhanced 
solution fidelity when compared to the Beam-Warming methods. The 
flux vector split technique as developed here follows the approach of 
Reklis and Thomas (Ref. 15). 

This technique, as with the Beam-Warming approaches, is based on 
a local linearization of the governing flow equations, Eq. (18). 
Equation (18) is locally linearized by employing the Jacobian matrix, 

J = 3F/8U, and the homogeneous property of the flux vector "f. The 
Jacobian matrix 3F/9U is diagonalized by a similarity transformation 
(Ref. 14) as 


where the matrices S and J are defined in Ref. 14. Decomposition of 
the Jacobian matrix is achieved by 

|£ = (u-a)4' + uf + (U+a)4' 

9u 1 2 ' ' 3 (45) 


u-a 


u+a 


(44) 


where the values are the decomposition functions (see Ref. 14). 
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Applying Eq. (48) to the governii.g Eq. (22), the difference equation 
written in flux split notation becomes 


[l + At (|u [Fi(a + |£) + Fa{a + |^) + Fa(a + |^)] 


-3 V 


3U 


Fx(a - |i) + Fa(a - |£) + ^3(3 _ Aijj _ gjj 


AU 


(49) 


= -At |Fi(a + 1^) + F2(a + |^) + F3(a + |^)] 


-[Fi(a - 1^) + p2(a - 1^) + Fa(a - |^)] - 6 - G^At| 


In this notation, the Jacobian matrices contained in Eq. (49) are opera- 
tors on the vectors AU(a), AU(a + A?), and AU(a - A?), thus providing 
the block tri-diagonal structure of the solution scheme. 

^ The method used to form the left hand side coefficients A, B, and 
C for the governing equation, which is in the same form as Eq. (26), 
is based on the extrapolation from cell faces to nodes based on the sign 
of the eigenvalues. The following sketch serves to illustrate this 
extrapolation process. 
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The eigenvalues, X*** and X-, as used here, denote positive and negative 
values, t'espectively . The sign of the eigenvalues, X, determines, for 


the flux piece at the cell face a - or a + A ^ , whether the 


O 
— L. 


extrapolation fills the A or B matrix or B or C matrix, respectively. 
For the case where all three eigenvalues are of the same sign (positive, 
for example), then matrices A and B are filled with the C contribution 
set to zero. If an eigenvalue changes signs within the three points. 
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i * h h i - "\, then the B component is set to zero, and the A and C 
matrices are filled. The Jacobian pieces required in Eq. (49) are 
evaiuated from the flux piece definitions given in Eq. (46). 

Boundary conditions for the split flux scheme are treated in a 
manner similar to the Beam-Warming boundary conditions except that 
Eqs. (29) and (30) are altered to characteristic form. Following the 
diagonalization of the Jacobian matrix as given in Eq. (44), Eqs. (29) 
and (30) are rewritten as 


aUj_j + SJB aUj * SJr* 

(50) 

SJB AUi + Sjr AU 2 = SJrt 

(51) 


where each equation in the above system is associated with an eigen- 
value (u - a, u, u + a) . The characteristic directions are used to 
determine the equation to be replaced with algebraic constraints. 

3.4.4 Split Characteristics Algorithm 

The Beam-Warming and split flux methods presented in the pre- 
vious subsections utilized implicit finite-difference techniques to solve 
the governing flow equations. Both approaches resulted in a block 
tri-diagonal solution structure requiring the inversion at each time step 
of a system of block tri-diagonal equations to yield the flowfield solution 
at the node points. These impiicit solvers, as shown in later sections, 
yield accurate solutions for complex flows including variable source 
terms over a “wide range of Courant numbers. 

The implicit algorithms have two basic problems, however, when 
applied to inlet flowfields. The first problem area lies in the basic 
solution technique which relies on the temporal linearization of the flux 
vector and source terms. This results in a relatively low maximum 
Courant number limit for complex flows. The second problem area lies 
in the execution time requirement of these implicit methods. Since a 
block tri-diagonal matrix must be inverted at each time step, solution 
times can become large for grids with a large number of node points. 



To partially overcome these two drawbacks of the Beam-Warming 
and split flux algorithms, split characteristics and MacCormack algo- 
rithms are presented. These two techniques are based on a predictor- 
corrector solution approach (Refs. 16 and 17) that yields the same 
spatial accuracy as the Beam-Warming and split flux schemes, but does 
not require block tri-diagonal matrix inversion. This results in faster 
overall solution times for the predictor-corrector schemes when com- 
pared to the block tri-diagonal based algorithms. In addition, it will be 
shown that the split characteristics technique can be extended to higher 
Courant numbers (resulting in additional savings in computer solution 
time) than the split flux or Beam-Warming methods for a given complex 
flowfield. 

The basic philosophy of the split characteristics approach is to use 
characteristic theory to determine the direction and magnitude of infor- 
mation propagation, which is accomplished by calculating the character- 
istic directions of the equations. The compatibility relationships, 
however, are not used in the conventional total derivative form. They 
are retained and used in partial derivative form. A conservative form 
of the governing equations is used to calculate the solution delta, AU. 
Characteristic theory is employed to split these solution deltas which 
are then propagated in the appropriate direction and transformed back 
to original variables. 

The result is a conservative, shock-capturing method using char- 
acteristic information. Boundary points are treated in a similar manner 
by combining the appropriate split components with known boundary 
conditions. This eliminates the need to arbitrarily guess which equa- 
tions to keep and combine with the boundary conditions. Character- 
istics theory determines the information to be retained and discarded. 

The following subsections contain a description of the split char- 
acteristics and MacCormack algorithms as applied to the quasi-one- 
dimensional flow equations. The subsections are devoted to the 
definition of: 

1 . Characteristic relations 

2. Explicit split characteristics algorithms 

3. Boundary conditions 
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4. Time integration algorithm 

5. Implicit operator 

6. Algorithm application 

3.4.4. 1 Characteristic Relations 
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This subsection is devoted to the development of the characteristic 
equations. Compatibility equations and characteristic directions are 
defined corresponding to the eigenvalues of the system of equations. 
These relations are required in order to accurately and clearly present 
the split characteristics concept. 

The governing system of equations defined in Eq. (18) can be 
written as 





- •» 

r.‘ 






U + * G* 

T X ^ 


(52) 


where G* = G - (5 )>.F and A = 3F/3U. The characteristics are obtained 

X ^ 

by diagonalizing A (represented in Eq. 44). As used here, the diagon- 
alizing matrix is given in Eq. (53). 


(53) 


The second matrix, J represents a transformation from conservative 

cp 

variables U to primitive variables Up where 




' 0 

-pa 

1 


1 

0 

0 


s = 


0 

-1 


-u/p 

1/p 

0 

i - v' * 

r'r, -.r . , 

— V‘^ 


0 

pa 

1 


(y-1)uV2 

-(y-1)u 

Y-1 


U = 


P 

pu 

E 


and 


p 

u 

P 


(54) 



The first matrix, then, transforms from primitive to characteristic form. 
This factored form is a convenient and efficient way to calculate the 
transformation. It is also needed for boundary condition calculations. 
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By diagonalizing the Jacobian A using Eq. (44), the governing 
equations can be cast into the partial differential form of the com- 
patibility equations 

■ ^ 1 ®* ( 55 ) 

where the compatibility equations hold along characteristic directions 
defined by 

^ (56) 


Equations (55) and (56) are the compatibility equation and char- 
acteristic direction pair corresponding to the eigenvalues X.. A typical 
characteristic method would finite difference the three compatibility 
equations along the corresponding three characteristic directions to 
calculate the U vector at a new time step. This requires interpolation 
since the characteristic directions do not necessarily align with the 

grid. Shocks also require special treatment in a characteristics method 
approach . 

3. 4. 4. 2 Explicit Split Characteristics 

It is the goal of split characteristics to use the characteristic 
information in a conservative form without requiring interpolation. In 
the split characteristics formulation, the partial differential form of the 
compatibility equation, Eq. (55), is used. The characteristic directions 
are applied, not as a direction for calculating total derivatives, but in 
determining the directions to calculate one sided partial derivatives. 
The following subsection separates the compatibility equations into 
pieces according to the characteristic directions, then transforms the 
resulting pieces back into conservative forms. 

In ord^er to split Eq. (55), it is convenient to define diagonal 
matrices D and D where the diagonal components, d, are 

d| = 1. 
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and D = I - D . Multiplying Eq. (55) by D“, Eq. (58) is obtained 


D"SU + » D"SG* 

X X ^ 


( 58 ) 


where A = SAS'^ with A and S defined in Eq. (52) and (53), respective- 
ly. This operation splits the system into two systems. Those obtained 
from D are ‘■he partial differential form of the non-negative 
characteristic values, while those obtained from D" correspond to the 
negative characteristic values. Each system contains a subset of the 
original three equations, and each equation is in one system or the 
other. 

By defining the matrix 1" = S"^D*S and noting that S-^D*AS = l‘A, 
the compatibility equation may be rewritten in the form given below. 

C59) 

Based on the form of Eq. (59), it is apparent that the operator 1" are 
splitting operators which split the original conservative equations into 
two parts, depending upon the characteristic directions. Equation (59) 
may be combined and rewritten as 

• (l‘ t nu^ = l‘tG - t rCG - (5/)^] (60) 

The basic principle behind the split characteristic method is that 
for the I operator, G - (^^F)e is evaluated using backward deriva- 

tives, while for I forward derivatives are used. Actually, G - (^ F)j. 
is evaluated or'y once per interval, as it is both a forward dii.erence 
for the right point and a backward di'Fference for the left point. 

From Eq. (60) for each node interval (at a cell boundary), an 
unsplit solution delta (AU) is applied as follows 

AU = At[G - (61) 


$ 
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where A? = 1 and j is the node number. To retain accurate wave speed 
characteristics, the pressure area term (pA^/A) in the source vector G 
is evaluated as 


pAj/A . (Pj^, t pj) (Aj^j - Aj)/(Aj,i + Aj) 


All other source terms are evaluated using averages of the interval end 
points. 

Once the unsplit deltas, following Eq. (61), are evaluated, the 
deltas are then split using the 1" operator defined in the previous 
paragraphs. The 1' operators are evaluated using a transformed 
average following Roe (Ref. 18). Splitting is performed per the fol- 
lowing matrix operation 

au* = i"AU = s"'d-sau 



The split AU" are combined at each point by summing the AU^ from the 
interval left of the point and AU from the interval right of the point. 

Time integration of the governing equations is based on a two- 
step, second-order Runge-Kutta algorithm. The method applied here is 
identical to the MacCormack explicit algorithm (Ref. 16) except for the 
manner in which the spacial averages are evaluated. Application of the 
two step Runge-Kutta algorithm gives 

U = + au 

= 1/2(U^ + U + aU) 

where represents the old or time level N known solution, and AU 
represents the combined split deltas obtained from Eq. (63). The 
superscript tilda designates predictor values. Time-dependent source 
terms are evaluated at the old time level and held fixed for the new 
level calculation. This method of treating the source terms is com- 
patible with that used in the split flux and Beam-Warming schemes. 




Application of Eq. (64) directly, with no splitting of the deltas, 
results in the explicit MacCormack algorithm. Fourth-order dissipation 
is employed for application of the explicit MacCormack method to eli- 
minate the formation of expansion shocks. 

3. 4. 4. 3 Boundary Conditions 

For application to supersonic inlets, the inflow boundary is 
assumed to be supersonic at all times. For supersonic flow at the 
inflow boundary plane, all three eigenvalues are positive. As a result, 
all of the characteristic directions point downstream. No information is 
passed upstream from the flow to the inflow boundary. Thus, all inflow 
conditions are specified from external considerations and are not part of 
the internal algorithms. As with the source terms, inflow boundary 
flow variables are needed at the old time-step and are held fixed for 
both the predictor and corrector steps. 

The flow at the outflow boundary is assumed to be subsonic and 
outflowing. Under these conditions, two characteristics (u*a and u) 
point downstream, while the (u-a) characteristic points upstream. The 
u-a characteristic information must therefore be discarded and replaced 
with a boundary condition constraint. As with the Beam-Warming and 
split flux solution schemes, boundary condition constraints implemented 
in the split characteristics scheme are specified static pressure, mass 
flow, and Mach number. A boundary condition is combined with the 
information from the u+a and u characteristics to obtain a AU for the 
outflow boundary point. 

For the outflow boundary, positive characteristic components are 
given by 



0 


0 

0 

0 ^ 


1 

0 

0 

P2 

s 

a2 

0 

-1 


-u/p 

1/p 

0 

Q3 


0 

pa 

1 


(y‘1)uV2 

-(y-1)u 

Y-1 

..a 


( 65 ) 
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where the characteristic components of AU, namely Q, are calcu- 
lated by Q = S AU. The above equation is used to calculate Qj and Q 
from the unsplit AU calculated from the last interval. The trans- 
formation matrices are evaluated using transformed averages of the last 
two points. 

The transformation from AU to Q, and is singular. That is, Q, 
and Q 3 may be calculated from AU, but AU cannot be uniquely calcu- 
lated from Qj and Q 3 . The additionally supplied boundary condition 

must be used to complete the system of equations and give a unique 
AU. 

Written in terms of the primitive variables U defined in Eq. ( 54 ) 
Eq. (65) becomes ^ 

0 0 0 Ap 0 

0-1 AU = Q 2 (66) 

0 pa 1 Ap Qa 

The system of equations can be solved for AU^ once the first equation 

in Eq. ( 66 ) is replaced with the outflow boundary condition or con- 
straint equation. 

The specified pressure boundary condition can be written as 

Ap - Ap* where Ap* is the known pressure increment. Replacing this 

eauation with the first equation in Eq. ( 66 ) and solving for AU results 
in P 


Ap * (Q 2 + Ap*)/a^ 

AU * (Qa + Ap*)/pa 
Ap = Ap* 

(67) 

The primitive variables 
by 

are related to the conservative variable. 

AU » AU 

cp p 

( 68 ) 


where is the inverse of as given in Eq. ( 53 ), 
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The specified mass flow and Mach number constraints are imposed 
in a manner similar to the pressure constraint except that linearized 
constraints equations are employed, For the mass flow boundary condi- 
tion, the linearized constraint equation is 

UAp + pAU = A(pu)* (69) 

with A(pu)* designating the known mass flow constraints. Utilizing 
Eq. (69) as the first equation in Eq. (66), the primitive variable con- 
straints become 

Au = (a2A(pu)*/0 - Q2-Q3)/(^aV0 - pa)^ 

AP * Qs - paAu 
Ap » (Qa + Ap)/a^ 

The ^ and^ variables are evaluated at the outflow point and not 
averaged as with the other elements of the transformation matrix. 

For the specified Mach number constraint, the isentropic Mach 
number relation is used where 


= puVp (71) 

Equation (71) is linearized and substituted in Eq. (66) leading to the 
primitive variable constraints given in Eq. (72) 

AP = [y$(aM^)* - - 2^0i^^Q3/pa]/ 

- 2^2^ Vpa - Y$0^) 

(72) 

Ap = (Qa + Ap)/a* 

AU = (Qa - Ap)/pa 

The conservative form of the variables is obtained by applying Eq. (68) 
as before. 


3. 4. 4. 4 Implicit Operator 

The implicit operator used in the split characteristic scheme is an 
artifice used to raise the Courant number limit of the explicit algorithm 
while maintaining the original order of accuracy. It can be thought of 
as a smoother of the explicit calculations. It redistributes the explicit 
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deltas using the convective behavior of the original quasi-one- 
dimensional equations. It is equivalent to an implicit method where only 
selected terms are evaluated implicitly. This partial implicit concept 
results in considerable computational savings over a fully implicit 
method. 

The operator used here is similar to MacCormack (Ref. 17) except 
that modifications are made to include split characteristic considerations. 
These modifications are easily made since the MacCormack algorithm 
already uses characteristic logic. 

Following the MacCormack concept, the governiiig equations, 

Eq. (18), are differentiated in time leading to 

“tr * (Sxf'st ■ (73) 

Reversing the order of differentiation and expanding the second term 
for a stationary grid gives 


U 

XT 




G 

T 


(74) 


Following the work of White (Ref. 19) the term in Eq. (74) is dis- 
carded, which is equivalent to treating the source team explicitly. The 
inclusion of these source terms in the implicit algorithm would require 
knowledge of the partial derivative of the terms with respect to the 
dependent variables resulting in an increase in the computational effort. 

Applying temporal differencing to the first term gives 

(u; - uj) * • 0 


where o and 0 represent old and new time levels, respectively. 
Splitting the Jacobian matrix. A, results in 


uj t 4t[(5/u8)j+ ■ U“ 


(76) 
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Using spacial differences of 5 over their domain of influence and 
replacing U® and with finite difference forms, Eq. (76) can be 
rev/ritten as 

+ [I + At 5 ^(A^ - A")] 6 Uj + = AUj ( 77 ) 

cfc B 

where U“ = AU/Ai and = 6 U/At has been used. Here j designates 

the centered node point and and A“ are considered operators on 6U. 
Equation (77) can be cast in a more symmetric form as given below 

-AT5x!A‘^l5Uj.l + U +At^|Al) 6 Uj - At5x1A“|6Uj+1 = aU^ ( 73 ) 

which is, in form, the same as the tri-diagonal system developed for 
the implicit Beam-Warming method. 

To further simplify algorithm solution requirements, the tri- 
diagonal implicit operator is factored into two bi-diagonal operators. 
The bi-diagonal systems, in addition to being simpler than a tri- 
diagonal system, easily convert to scalar equations employing the 
characteristic transformation used for splitting. The factorization is 
written as 

(I + AtcJA^DaUj = AUj + ATCjA‘^|6Llj_i (79) 

(I + AtCj^IA"! ) 6 Uj = 6 Uj + AtC^IA (80) 

The factorization is considered approximate as a higher order (At)* 
term is dropped . 

The characteristic transformation used to split the implicit operator 
is used to simplify the solution of the bi-diagonal systems. As applied 
here, these operators are actually block bi-diagonal systems with 3x3 
blocks. The characteristic transformation provides for a solution to 
this system. 

It is convenient to describe the algorithm assuming the right hand 
side (RHS) of the operators is already calculated. The algorithm then 
has two steps: 


36 






1. Solve for the and 6U. following Eq. (79) and (80). 

2. Propagate the appropriate terms to the RHS of the next 
spacial point. 

One of the advantages of the MacCormack algorithm (Ref. 17) is 
its ability to be implicit where needed and to be explicit otherwise. 
This is accomplished by shifting the diagonal values of the implicit 
operator. The same principle can be used with the split characteristics 
scheme presented here. Denoting |X."| as the diagonal values of |A‘|, 
the parameter X.* can be defined as 

A+ X 

Xj- = max (0. |x^| - CFAj^At) (81) 

The split characteristics algorithm uses the magnitude and sign of X for 
defining |A | instead of the magnitude of the eigenvalue as used in the 
MacCormack method. 

The Courant Factor (CF) is the Courant number level at which the 
implicit operator is applied. Since the present explicit split algorithm 
appears stable to Courant number values slightly above 1, the CF was 
set at 1 . The result is an algorithm that is explicit at points with a 
local Courant number less than 1, while smoothly transitioning to an 
implicit node at points where the Courant number is greater than 1. 

3. 4. 4. 5 Algorithm Application 

With the typical tri-diagonal methods, the boundary conditions and 
the tri-diagonal operator are combined into one system to be solved 
together. For the split characteristics scheme presented here, the 
bi-diagonal systems and the boundary conditions can be completely 
separated. However, there is an implied order of calculation. It is 
important to emphasize here, that a supersonic upstream boundary is a 
critical part of the factored algorithms. This supersonic condition 
separates the upstream boundary condition from the flow. Character- 
istic theory shows that the internal flow cannot influence the upstream 
boundary point. It is this fact that breaks the implicit coupling of the 
implicit operators and allows a starting point. 
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Application of the factored implicit algorithm can be described in 
the following four steps. 

1. The explicit split deltas are calculated following Eq. (63). 

2. The first implicit operator as given in Eq. (79) is applied by- 
starting with the first interior point and ending with the 
outflow boundary point. This operator passes information 
downstream and is called the forward implicit operator. Since 
tlje upstream boundary condition is externally determined, the 
6lJ can be calculated in order, each point using values -from 
the previous point. 

3. The forward implicit operator also operates upon the partial 
AU calculated at the outflow boundary point. The outflow 
boundary roy^^ne is identical to the external routine except 
that it uses 6U instead of AU. 

4. After applying the appropriate outflow boundary condition, 
the backward operation as given in Eq. (80) is used. This 
operator starts at the last interior point and ends with the 
first interior point. Since |B~| is null where flow is super- 
sonic, the operator is a trivial identify operator for the 
leading supersonic points. The results of the outflow 
boundary algorithm are used to start the backward operator. 

5. The variables are updated with the same second order Runge 
Kutta algorithm using 6U instead of AU as implied in 
Eq. (64). 


3.5 GEOMETRY 

Capabilities are provided for three types of axisymmetric or two- 
dimensional inlet geometries. The first geometry type is a simple duct 
with no centerbody. This geometry is useful primarily for comparison 
with other numerical solutions. 

The second geometry type adds a fixed contour translating center- 
body which may have one or two straight segments upstream of the cowl 
lip in the axisymmetric case or up to three in the two-dimensional case. 
The centerbody translates forward from a specified point to provide 
restart capability. A sketch illustrating this geometry is shown in 
Fig. 3. 

The third geometry type has a variable diameter (or height) 
centerbody with a slot for bleed. Two straight segments upstream of 
the cowl lip are specified. The second segment is allowed to effectively 
rotate about the point joining the first two segments. A slot follows 
the second segment. The section downstream of the slot is also allowed 
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FIGURE 3. FIXED CONTOUR CENTERBODY GEOMETRY 




FIGURE 4. VARIABLE DIAMETER CENTERBODY GEOMETRY 
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to effectively rotate. In addition to the primary centerbody translation, 
a secondary translation of the first two segments is provided. A 
shetch illustrating this geometry is shown in Fig. 4. 

The inlet cowl, centerbody, and in the two-dimensional case, width 
contours, are defined by up to 20 segments. The segments may be 
linear, parabolic, or cubic except for the centerbody segments upstream 
of the cowl lip which must be linear. The segments are defined by 
specifying end points and, in the case of parabolic or cubic segments, 
intermediate points or slopes at the end points. 

3.6 UNSTART MODEL 

When a supersonic inlet unstarts, the normal shock wave is 
expelled from the inlet and moves upstream toward the vertex of the 
centerbody as illustrated schematically in Fig. 5. Behind the normal 
shock wave the flow is subsonic, and, because the normal shock wave 
is detached from the inlet lip, there is spillage of the incoming air over 
the inlet housing. The operation of the inlet diffuser under this con- 
dition is said to be subcritical. 



FIGURE 6. INLET UNSTART SCHEMATIC 
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Modeling of the above discussed inlet unstart pheno-<...na is herein 
accomplished through adaptation of an incremental control volume 
approach. The Moeckel continuity method (Ref. 3) Is used to relate the 
expelled normal shock position ahead of the inlet cowl to the amount of 

spilled mass flow over the Inlet housing. Details of this approach are 
explained below. 

Appendix A presents a methodology whereby the actual spilled 
mass flux can be determined in terms of the straight sonic line mass 
flux and a correction factor F which is a function of free-stream Mach 
number only following Moeckel (Ref. 3). This approach is implemented 
as Illustrated in Fig. 6 where the incremental control volumes are 
formed by the x-direction one-dimensional numerical grid mesh. 
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The actual spilled mass flux pV from each local control volume is 
determined via a modified form of Eq. (A-3) In Appendix A as 


i(- , 


; ; V. •'•'‘.'’.(L' , 

•-'VMviA' ' 




pV X p*v» 


(82) 


with the sonic mass flux p*V* given by an isentropic expansion from 
the local control volume stagnation (total) pressure and stagnation 
(total) temperature to the local sonlc(*) conditions. The correction 
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factor is a modified form of the straight sonic line correction factor 
F derived in Appendix A and accounts for inlet unstart variable area 


= F 
c 


Fg is given by 

- 



•’A / \\J 



Wc ; [-ij 

1 - 



PlUl ^ 

4 


f(83) 


V- Ui. 

where j=0 for two-dimensional geometries and j=1 for axisymmetric 
geometries. The geometry parameters y^, a^^, and 


3j are defined In 


Fig. 7 where it is noted that for the special condition y^ = aj/a^^ 
Eq. (83) reduces to F^=F which is consistent with the derivation of the 
straight sonic line correction factor F per Appendix A. 



PIQURB 7. MODIFIED CORRBCTION FACtOR QEOMETRV 
PARAMETER DEFINITION 

Numerical implementation of the above discussed incremental control 
volume approach is accomplished by appendage of additional grid points 
in front of the inlet cowl upon expulsion of the normal shock. The 
moving normal shock location is determined by monitoring the computed 
Mach number distribution and interpolating for the Mach one location in 
the shock jump from supersonic to subsonic. Note that in this one- 
dimensional approximation, the moving normal shock is assumed normal 
to the centerbody x-axis coordinate. Determination of the straight 
sonic line correction factor F foi use in Eq. (83) is per Fig. (A-1) in 
Appendix A, with the free-stream Mach number taken to be the local 
average Mach number in the centerbody inviscid flow field immediately 
upstream of the normal shock, i.e., in Fig. 7. The mass flux pV in 

Fig. 7 is determined in the exact same manner. F and F are assumed 

o 

constant for all control volumes at the values determined above. 
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Also considered in the incremental control volume analysis are both 
the momentum flux and energy flux spillage in a manner exactly ana- 
logous to the mass flux spillage presented above. Even though the 
momentum flux is, in reality, transferred out of the system in a normal 
direction, one-dimensional gas-dynamic theory requires that the momen- 
tum be accounted for as if the transfer were accomplished in the tan- 
gential direction. This is necessary in order to satisfy the Second Law 
of Thermodynamics and is discussed on page 234 of Shapiro’s book 
(Ref. 8). Basically, normal mass extraction generates an entropy 

decrease due to an ambiguity in the quasi-one-dimensional inviscid flow 
equations. 


With ihe mass, momentum, and energy flux spillage for each 
individual control volume determined as discussed above, unstart 
phenomena are modeled through the governing equations of motion (1-3] 
source/sink terms M^, F^, and Q^. For example, mass spillage due to 
inlet unstart is treated locally as a mass sink and incorporated into the 
continuity equation (1) source/sink term M^. In this manner all of the 
fluid dynamic effects of inlet unstart are properly incorporated into the 
numerical solution procedure in a totally consistent manner which 
reflects the local conditions of the flow field. It should be noted that 
this quasi-steady approach also properly treats combined phenomena 
which may occur locally, such as mass loss due to spillage and mass 
loss due to centerbody bleed in the same control volume. In this case 
the continuity equation (1) source/sink term would be the sum of 
the separate spillage and bleed within the control volume. Similar 

comments apply to the momentum and energy equation source/sink terms 
Fg and Q^, respectively. 


3.7 BLEED AND BYPASS MODEL 

Inlet bleed and bypass flow rates are computed using the local 
stagnation (total) pressure and stagnation (total) temperature of the 
inlet duct flow in conjunction with bleed and bypass flow coefficients 
from experimental data. For cases where the bleed or bypass flow is 
not choked, the bleed or bypass plenum conditions (pressure and 
temperature) are required to determine the flow coefficients. These 
unchoked plenum conditions a»*e computed by assuming the plenum to be 
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a single lumped volume with negligible flow velocity and numerically 
solving two ordinary differential equations for the time-dependent 
plenum density and temperature; plenum pressure is then determined 
through the equation of state. The plenum discharge conditions are 
assumed to be choked flow. A schematic of an inlet bleed system is 
shown in Fig. 8; an exactly analogous situation applies for an inlet 
bypass system with bleed flow replaced by bypass flow. 
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FIGURE 8. INLET BLEED SYSTEM SCHEMATIC 
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The governing differential equations for the bleed and bypass 
model are taken from Chapter 2 in the textbook by Zucrow and Hoffman 
(Ref. 9). Assuming uniform properties inside the constant volume 
plenum, continuity of mass requires 


dptj 1 

* i(F~ “ PzUzAa) 

while conservation of energy requires 


(84) 
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Under the assumption of negligible flow velocity within the plenum 

(Up = 0) as well as conservation of total enthalpy across flow 
boundaries, i.e., 

' Cp<''T>0uct 

Eq. (85) can be written as 


dt V p 
P^P 


■ (Y-l)TpP2U2A2 


With the substitution of Eq. (84). Given the plenum density and 
temperature determined via numerical integration (trapezoidal rule) of 

Eqs. (84) and (87), the plenum pressure is calculated from the equation 
of state (4) as 

P_ * P_RT_ 

P P P (88) 

Experimental data from Dennard (Ref. 20) and Syberg and Hickcox 
(Ref. 21) are used to determine the bleed flow coefficients as a function 
of bleed hole geometry parameters and local inlet duct flow conditions. 
Unpublished flow coefficient data are used for the bypass flow coef- 
ficients. The bleed or bypass flow is determined from the flow coeffi- 
cient relation 

Pt Ai 

PiUiAi = Ki 7=i= Fun(v) 

where it is assumed that 

Pt, “ ^Pl^Duct 

ll I UUCt 

hx ’ ^^T^Duct 

and Fun (y) is a function of the specific heat ratio only. 





The flow coefficient Ki is the ratio of the actual mass flow through the 
porous bleed section (or bypass) to the maximum theoretical value for 
choked flow. It is a function of the following flow and geometry 
variables: 

• Hole length/diameter 

• Boundary-layer thickness to hole diameter 

• Hole angle 

• Hole area to total area (porosity) 

• Local Mach number 

• Pressure ratio across hole 

The present work, limited to inviscid flow considerations, utilizes a 
table of values for Ki as a function of hole geometry with local Mach 
number and pressure ratio as parameters. For the choked flow plenum 
discharge, the mass flow is given by 

p? A2 

P2U2A2 * K2 Fun(Y) (9D 

where it is assumed that 


■ Tp C92) 

where recall that the flow velocity within the plenum volume is assumed 
negligible. The flow coefficient Kj is the choked mass flow coefficient 
for the plenum exhaust valve which has a typical constant value 
between 0.9 and 1.0. For the present work, K* is taken as constant at 
the value 0.9. 

Initial conditions from which to start the numerical integration of 
Eqs. (84) and (87) are based upon the condition of steady-state flow 
through the plenum. Under this restriction, the governing equations 
become 

PlUlAx sa P2U2A2 ronx 


^’’’r^Duct * "’’Ti “ ’’’p “ ''^12 
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4. Using known conditions at the cone surface and the shock, 
approximate the local flow angle 0 distribution through a 
cubic fit 

# * t - cos"^ (y) 


u * ai + + a3(tj)-tj)^)® 


V- 3a 3 (<i< -<!<*) 

“ 2(V^ “ —2 

2V^ 2V„cosi)»5 

5. Compare the integrated mass flow between the cone surface 
and the shock with the corresponding free-stream capture 
mass f<ow. Iterate to convergence using the secant method to 
deter'rnine the updated shock angle 'j»g. 


(99) 

(ino) 

( 101 ) 

( 102 ) 

(103) 



FIOURE 9. INVI8CIO CONICAL FLOW riELO NOMBNCLATURC 

Figure 10 represents a comparison of the present approximate 
technique relative to exact results from the inviscid conical method, 
Zucrow and Hoffman (Ref. 9). Here the flow field is for a 10.0 degree 
semi-vertex angle sharp cone in a Mach 2.50 free-stream flow which is 
representative of an actual mixed-compression inlet centerbody. Agree- 
ment between the approximate and exact techniques is excellent for all 
flow field variables and shock location. 
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For the case of a mixed compression axisymmetric inlet with a 
biconic (fore-aft cone) centerbody, the supersonic external inviscid flow 
field contains an embedded shock, and the above discussed approximate 
flow field technique must be modified to reflect this phenomena. With 
reference to Fig. 11, the following procedure is utilized in the 
present work to incorporate an embedded shock capability into the 
approximate conical flow technique: 

1. The initial embedded shock angle is determined from classical 

oblique shock theory based upon the fore cone surface Mach 

number and the compression angle *l >2 - > This shock angle 

c 

is used to determine the shock point "1." 

2. Based upon examination of numerous- inviscid method of char- 
acteristics solutions (using Ref. 26) for various biconic 
configurations in supersonic flow, it is assumed that the total 
pressure loss across the embedded oblique shock is the same 
(i.e., constant) for ail embedded shock points. This total 
pressure loss is determined via the initial shock angle dis- 
cussed above. 

3. With embedded shock point "1" known, determine embedded 
shock point "2" from the known (constant) total pressure loss 
and the local upstream flow conditions on the corresponding 
conical ray. 

4. Repeat this formal procedure to "march" the embedded shock 
points out to point "s" where the embedded shock intersects 
with the conical shock from the fore cone. 

5. Assume the inviscid flow field behind the embedded shock is 
conical with respect to the virtual origin of the aft cone. 

6. "Correct" the local flow field pressures behind the embedded 
shock by multiplying by the ratio of the aft cone surface 
pressure (based on inviscid conical flow theory) to the two- 
dimensional surface pressure (based on oblique shock 
theory) . 

7. Compute the remainder of the aft cone inviscid flow field 
variables using the assumed conical flow field condition in 
conjunction with an isentropic expansion from conditions 
immediately behind the embedded shock to the local 
"corrected" pressure on an aft-cone conical ray. 

Figures 12 and 13 show comparisons of the present appt'oximate 
technique relative to exact results from the inviscid method of char- 
acteristics program of Inouye, Rakich, and Lomax (Ref. 26). Here the 
flow field is for a biconic configuration with a 10.0 degree semi-vertex 
angle fore cone and a 18.5 degree aft cone in a Mach 2.50 free-stream 
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flow which is representative of an actual mixed-compression 
inlet conterbody. As can be seen from Fig. 12, location of the 
aft-cone-generated embedded shock is reasonably well determined by the 
approximate technique with the largest discrepancy near the intersection 
with the fore-cone shock. The comparison of flow field profiles 
(velocity, density, and static pressure) across the embedded shock 
layer at x-location of 2.12 is given in Fig. 13, where 10 conical ray 
intervals have been used across the aft-cone shock layer. Agreement 
between the approximate technique and the exact method of character- 
istics solution is good across the entire shock layer, both in distri- 
bution and magnitude. 

With the profile distribution of flow field properties (vel <city, 
density, and static pressure) across the inviscid shock layer available 
from the above discussed approximate techniques, it remains to deter- 
mine equivalent one-dimensional quantities for use as upstream boundary 
conditions at any plane (see Fig. 14). 



FIGURE 14. UPSTREAM BOUNDARY CONDITIONS 


Plane "A" is the appropriate upstream boundary condition for a started 
inlet while plane "B" and "C" apply to unstarted inlets with an expelled 
normal shock following the unstart model presented in Section 3.6. 
Formal integration across any plane normal to the x-axis yields 

Qi = y pudA 
A 


Mass 


(104) 
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p 2 » y (p+puV)dA 
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Momentum 


(105) 


(106) 


Qs * y (E+p)udA Energy (106) 

A 

where V is the total velocity along a streamline and u is the corre- 
sponding velocity component parallel to the x-axis. In terms of the 
basic flow variables Uj, Uj, and U, 


corre- 


Qi ~ ?vL|2 


(107) 


Qa = [(y+1)U3 + (-^) 


(108) 


3. . Cvu, - (0^) ^ 


(109) 


where Q = x^Q. Solving Eqs. (107), (108), and (109) for Uj, Uj, 
U 3 yields 


= 1 


( 110 ) 


. -B ± Vb^-4AC 
U 2 2A 


(111) 
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The choice of sign in the quadratic root of Eq. (Ill) must be nega- 
tive (-) in order for the resulting flow to be supersonic; a positive (+) 
sign results in subsonic flo-, . 

In the case of a two-dimensional multiple ramp external flow, 
calculation of the flow field is much simpler since straight shocks are 
generated with uniform flows behind them. The classical inviscid, 
two-dimensional, oblique shock relations are applied using the secant 
method to converge an iterative solution process. By using an initial 
guess to the solution based on small disturbance theory (Ref, 25), 
convergence is generally obtained in two to four iterations. 

The flow field over a downstream ramp is computed using the flow 
over the preceding ramp as an upstream condition. Up to three 
straight ramps upstream of the cowl lip are allowed in the formulation, 

3,9 INITIAL CONDITIONS 

In order to provide a consistent flow field from which transient 
solutions may evolve, initial conditions are required within the inlet and 
upstream of the inlet for unstarted flow models that are compatible with 
the governing one-dimensional equations presented in earlier sections 
and with the engineering flow models that describe the engineering 
characteristics of the inlet configuration, both in the unstarted and 
started mode. The initial conditions must include nomal shocks within 
the inlet and nomal shocks for the unstarted inlet and also include the 
double-shock case for the unstarted inlet. Since most flow fields are 
established based on downstream boundary condition constraints, the 
initial conditions must use this downstream or exit plane boundary 
condition information to define the flow field within the inlet for both 
the unstarted and started inlet cases. Boundary condition constraints 
as indicated earlier, may be qualified as either pressure, Mach number, 
or mass flow constraints. 

The initial conditions are established within the inlet based on the 
steady-state one-dimensional isentropic and normal shock compressible 
flow relations. Using various combinations of these equations in a 
direct or iterative fashion as required by the boundary condition con- 
straints provides for a number of capabilities. These capabilities 
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include the steady-state flow field definition for both started and 
unstarted inlets which are coupled with the external flow field cal- 
culations. The unstarted inlet model used in the steady-state initial 
condition definition is consistent with the unstart model discussed 
previousiy. Consistent with the governing equation and engineering 
model formulations presented in earlier sections, the initial condition 
routines are generalized to include both two-dimensional or axisymmetric 
configurations. The initial condition capabilities provide the range of 
steady-state boundary conditions that define the allowable input bounds 
for the exit plane boundary conditions and other initial conditions. 
This range of boundary conditions is an especially important feature of 
the initial condition routine because the maximum and minimum input 
bounds are not always obvious. The initial condition routine does not 
include bleed or bypass modeling. This is not a serious drawback, 
however, since a given inlet configuration can be started in a mode 
without bleed or bypass. Bleed or bypass can then be "turned on" and 
3 new steady-state solution can be evolved in time* 

As stated in the previous paragraph, both started and unstarted 
inlet flow fields can be defined using the initial conditions steady-state 
routines based on the one-dimensional isentropic and normal shock 
compressible flow relations. Inlets that can be considered are both 
variable or con. tant area in configuration. For the started inlet con- 
figuration, exit plane constraints such as exit plane pressure, Mach 
number, or mass flow rate can be imposed or the shock location within 
a started inlet duct can be defined. All conditions assume, of course, 
that the flow field at the inlet plane for the started configuration is 
supersonic. The only restriction on the started iniet configuration is 
that the shock location must be downstream of the throat. 

For the unstarted inlet configuration, the freestream flow field 
must be supersonic in form with exit plane inlet constraints based on 
exit plane pressure, Mach number, and mass flow rate. As an added 
capability, the unstarted inlet configuration can be also specified by the 




shock location upstream of the inlet plane for the single-shock case and 
the upstream shock location and the shock location downstream of the 
throat within the inlet for the double-shock case. The only restriction 
for the inlet flow field for the unstarted double shock inlet configura- 
tion is that subsonic flow must be present at the inlet plane. 

The flow field for started inlets with the shock located downstream 
of the throat is based on the approach given in the textbook by Hall 
(Ref. 27). In this approach, the exit pressure, P 2 , divided by the 
upstream stagnation pressure, p-j-i, is defined as the exit plane 
boundary condition constraint. The exit plane Mach number, M 2 , can 
be determined by application of the equation given below. 

Y+1 


£2. . JL 

^ A* M 2 U+lJ I ^ 2 '^2 ] 


(116) 


This equation is derived by eliminating the entropy between the 
equations for the stagnation pressure for the upstream flow field 
divided by the exit plane pressure as a function of the exit plane Mach 
number and the exit plane area, Aj, as a function of the sonic area, 
A*, and exit plane Mach number, M 2 . The sketch given below shows a 
schematic of the region 1 and region 2 points for reference. 




Once the exit plane Mach number has been defined, isentropic expan- 
sion is used to relate the exit plane pressure to the local stagnation 
pressure in region 2. The Mach number just ahead of the shock front 
in region 1 is then defined by relating the region 1 Mach number to the 
ratio of the stagnation pressure in region 2 divided by the stagnation 
pressure in region 1. Once the Mach number in region 1 is deter- 
mined, isentropic area expansion can be used to define the cross- 
sectional area at the shock location as a function of the sonic area. 
Using the shock location as an internal constraint, the shock location is 
Iterated using the method of steepest descent to satisfy the exit plane 

boundary conditions whether they be pressure, Mach number, or mass 
flow rate. 

Unstarted inlet initial conditions are handled in much the same way 
as the started inlet flow field conditions except that the inlet entrance 
Mach number for unstarted inlets is subsonic. The flow field ahead of 
the inlet entrance and behind the upstream nomal shock is modeled 
following the unstart modeling concepts presented in previous sections. 
The grid cells in the region upstream of the inlet entrance are modeled 
using conservation of mass, momentum, and energy applied in a control 
volume sense. This leads to relationships describing the pressure and 
temperature jump across a control volume as shown in the sketch below 
between nodes j-1 and j. 


ORIGINAL r* i 
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The mass flow rate crossing the control volume boundary denoted here 
as r is related to the mass flow rate crossing a fictitious boundary 
whose properties are at the local sonic conditions defined by the 
Moeckel relationship as given earlier. Equations (117) and (118) 
describe the temperature pressure relationships for this control volume 
model which includes a mass bleed across the boundary (denoted with 
subscript r) . 


“j-i 

1 m. , m 

Pj “ Pj-1 ^ " r 


(117) 

(118) 


where 3 



In a manner similar to the started inlet case, the upstream normal 
shock location is iterated using the method of steepest descent to 
satisfy an imposed exit boundary condition whether it be pressure, 
Mach number, or mass flow rate. The double shock initial condition is 
established by iterating the shock upstream of the inlet plane until the 
inlet throat Mach number is exactly unity (M = 1.0). A shock down- 
stream of the sonic throat is then allowed to form to satisfy the exit 
plane boundary condition constraints. 


3.10 COMPUTER CODE (PROGRAM LAPIN) 

A schematic of the computer code (Program LAPIN - Large Perturba- 
tion Inlet) is shown in Fig. 15. The program itself is written in 
FORTRAN IV for execution on a Digital Equipment Corporation 
VAX-11/780 and contains approximately 6,100 lines of source code. All 
input is through NAMELIST with an example for the NAS.A LeRC 40-60 
inlet geometry given in Fig. 16. Output from Program LAPIN for a 
Mach 2.50 free-stream flow and a corrected mass flow exit plane 
boundary condition is reproduced in Figures 17 through 19. The initial 
(steady-state) distribution for an unstarted inlet is shown in Fig. 17 
while Fig. 18 gives the transient distribution at a time of 0.0240 
seconds. Note in Fig. 18 the internal embedded shock at x-location 
between 4.8625 and 5.0052 as well as the expelled shock at a x-location 
between 1.8649 and 2.0090. These calculations were performed using 
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FIGURE 15. SCHEMATIC OF PROGRAM LAPIIM (LArge Perturbation IMIet) 
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FIGURE 17. PROGRAM LAPIN INITIAL CONDITIONS OUTPUT 
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the implicit split characteristics algorithm with 41 grid points within th 
inlet (2.0090 S X S 7.71GO) and 7 grid points on the unstarted center 
body (1.0000 S X S 2.0090). The corresponding bypass plenum sum- 
mary output for this case is given in Fig. 19. 
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OF * 4 0 ALGORITHM VERIFICATION 


In order to assess the relative accuracy and solution fidelity of 
each of the algorithms discussed in the previous sections, an analysis is 
presented for the unsteady normal shock wave in a constant area duct. 
The five algorithms examined are the Beam-Warming, Hybrid Beam- 
Warming, split flux, split characteristic, and the explicit MacCormack 
solution scheme. The moving normal shock in a constant area duct 
problem was chosen as the basis for algorithm verification because exact 
analytical closed-form solutions are available for shock speed, as well as 
flow properties downstream of the shock front. 

By following Chapter 8 in the textbook by Owczarek (Ref. 28), it 
can be shown that the speed of the moving shock wave relative to the 
gas into which it is moving is given by 


with 



(119) 


( 120 ) 


where subscript "1" denotes upstream conditions into which the shock is 
moving, and subscript "2" denotes downstream conditions behind the 
shock per the following sketch. 


STATE AHEAD OF THE SHOCK WAVE 


STATE BEHIHD THE SHOCK WAVE 
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Auxiliary relations include 
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with the shock speed given by 
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For a given ratio of downstream to upstream mass flux ( p u ), 
the above equations can be solved iteratively to yield the moving shock 
speed as well as all downstream shock parameters. This provides an 
excellent test case for any one-dimensional unsteady inviscid numerical 
technique in that the moving normal shock in a constant area duct can 
be generated by a step-jump reduction in mass flow rate. 

Presented in Table 1 are calculated flow parameters from the above 
equations for the case of a Mach 2.0 upstream condition with a 
20 percent reduction in mass flow. Shown enclosed in ( ) are the 

corresponding stationary normal shock values for a Mach 2.0 flow. As 
is obvious from the table, the moving shock wave causes substantial 
changes in shock-jump conditions, including an almost seven percent 
increase in total temperature. 

Presented in Figs. 20 and 21 are the computed results for a 
Courant number of 0.7 at a time of 0.0050 seconds after flow initiation 
with a step jump reduction in mass flow at time zero as a prescribed 
downstream boundary condition. The upstream static temperature was 
taken to be 300°K (540®R) which results in the moving shock being 
located at a distance of 0.7365 based upon the shock Mach number 
given by Eq. (129), This distance location should be used by the 
reader to access shock speed performance of the various algorithms. 
As can be seen by comparing the results for the five different methods 
as described in the numerical algorithms section given earlier, all 
methods give relatively good shock speed results based on the analytical 
shock location of 0.7365 as indicated in the previous paragraph. 
Referring to these figures, it is difficult to determine the exact shock 
location for each method because of shock smearing. 

Some general conclusions can be drawn from these figures, however, 
to highlight the important features of the methods. Comparing the 
smoothness of each of the calculations, the split characteristics solution 



TABLE 1. UNSTEADY NORMAL SHOCK WAVE IN 
A CONSTANT AREA DUCT 
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FIGURE 21, MACH lUUMBER-DISTAMCE HISTORY AT TIME=0.005 SEC 
FOR AN UNSTEADY SHOCK IN A CONSTANT AREA DUCT 



scheme is shown to provide the cleanest or highest fidelity numerical 
solution for this shock tube problem. This solution exhibits no numer- 
ical or solution noise in contrast to the other four methods which do 
exhibit some form of numerical noise. The split flux algorithm solution, 
when compared to the split characteristic technique, gives nominally 
good results with some overshoot in pressure just downstream of the 
shock location. The hybrid scheme as compared to the split flux 
scheme gives nominally the same result except that the overshoot down- 
stream of the shock is somewhat greater. It must be remembered that 
the hybrid scheme results include a conservative fourth-order dissipa- 
tion which must be adjusted based upon a local Courant number condi- 
tion. The Beam-Warming and MacCormack schemes display significant 
solution noise, both upstream and downstream of the shock location 
when compared to the split characteristic, split flux, or Beam-Warming 
hybrid schemes. The upstream undershoot in the pressure for both the 
Beam-Warming and MacCormack schemes generally result in algorithm 
failure for stronger shock calculations. 

Similar sets of calculations have been performed for a para- 
bolic 2-1-2 nozzle configuration with an inlet Mach number of 2 and 
Courant number of 0.8 with comparable results for the shock tracking 
capability of each of the schemes. Calculations have also been made for 
the 2-1-2 nozzle to assess Courant number effects on shock tracking 
and fidelity capabilities of each of the numerical schemes for a ham- 
mershock case. These studies show that the split flux and hybrid 
schemes give little solution deterioration for Courant numbers as high 
as 6. The split characteristics scheme could be pushed, however, to 
Courant number values as high as 10 for the hammershock case. For 
reference purposes, the hammershock was moving at approximately one 
grid cell per time step at a Courant number of 12. Thus, both the 
split flux and split characteristics schemes can be utilized with con- 
fidence under hammershock transient conditions as long as the Courant 
number constraint is imposed providing at least two time-step calcula- 
tions per grid cell movement of the shock. 


5.0 RESULTS AND DISCUSSION 


5.1 40-60 INLET 

The NASA LeRC 40-60 inlet is an axisymmetric mixed-compression 
inlet with 40 percent of the effective supersonic area contraction 
occurring externally and 60 percent supersonic area contraction 
occurring internally at the design free-stream Mach number of 2.50. 
Details of the 40-60 inlet configuration and dynamic pressure transducer 
locations are given in Fig. 22. The inlet was equipped with a trans- 
lating centerbody, high-response sliding plate overboard bypass doors, 
and an ejector bypass used for engine cooling. Pertinent details con- 
cerning experimental investigations of this inlet conducted in the NASA 
LeRC 10- by 10-foot supersonic wind tunnel are contained in Refs. 29 
through 34. 

5.2 40-60 INLET STEADY-STATE PERFORMANCE 

The present section will computationally examine the steady-state 
operation of the 40-60 axisymmetric mixed-compression inlet at a free- 
stream Mach number of 2.50. Experimental data reported in Ref. 31 
will be used for comparison with calculated results. Details of the 
40-60 inlet configuration and dynamic pressure transducer locations 
have been previously given in Fig. 22. 

Overall inlet performance in terms of total pressure recovery as a 
function of engine-corrected mass flow variation is shown in Fig. 23. 
The calculated total pressure recovery is in good agreement with the 
experimental data and even predicts a flow oscillation at low mass flow 
which was experimentally observed at approximately the same value of 
mass flow. It is to be noted that the Program LAPIN momentum source 
term is based upon bleed and bypass extraction of momentum as if the 
transfer were accomplished in the tangential direction, even though the 
momentum flux is, in reality, transferred out of the inlet flow in a 
normal direction (see the discussion of this point in Section 3.6). If 
the momentum source term is not treated in this manner, erroneously 
high-pressure recovery (in some cases larger than unity) is predicted 
by the present approach. 
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FIGURE 22. 40-60 INLET DETAILS 
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Inlet static pressure distributions shown in Fig. 24 for the cowl 
surface and centerbody show that the calculated terminal normal shock 
location is significantly downstream of the experimentally indicated 
location. This discrepancy can be traced to the inviscid treatment of 
mass bleed effects without any allowance for viscous (boundary layer) 
interaction. Downstream of the throat region, the terminal normal 
shock and the boundary layer interact with possible flow separation and 
reattachment. Viscous-inviscid interaction is obviously dominant in 
establishing the resulting location of the terminal shock. Irrespective 
of this, the calculated static pressure distribution is in excellent 
agreement with experiment in regions upstream of the bleed locations 
and downstream of the calculated terminal shock. 

5.3 40-60 INLET HAMMERSHOCK TRANSIENT 

Hammershock transients are one of the most potentially dangerous 
situations which can occur in a supersonic inlet due to the large over- 
pressure levels involved. The present section will computationally 
examine a hammershock transient in the 40-60 axisymmetric mixed- 
compression inlet at a free-stream Mach number of 2.50. Experimental 
data reported in Ref. 33 will be used for comparison with calculated 
results. Details of the 40-60 inlet configuration have been previously 
given in Fig. 22. 

The experimental investigation of Ref. 33 generated a hammershock 
transient via engine compressor stall which resulted in a strong pres- 
sure disturbance, propagated upstream through the inlet. For pur- 
poses of performing the Program LAPIN calculation, the experimentally 
determined exit plane static pressure distribution was input as a pre- 
scribed downstream boundary condition. The results of this calculation 
are shown in Fig. 25, where station "2" represents the exit plane 
location. Inlet bypass setting was determined in order to match the 
experimentally determined total pressure recovery. The calculated 
upstream results are in reasonable agreement with the experimental data 
considering the fact that the present inviscid model positioned the 
steady-state terminal normal shock too far downstream due to neglect of 
viscous effects as noted earlier. 
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FIGURE 24. STATIC PRESSURE DISTRIBUTIONS 








It is of interest to computationally follow the development and 
movement of the above discussed hammershock transient through the 
inlet. Figure 26 gives the cross-sectional area distribution of the 40-60 
inlet where it should be noted that the throat is located at a x-value of 
approximately 3.5. Figure 27 shows the corrected mass flow distribu- 
tion and corresponding mass flow disti ibution; associated flow-field dis- 
tributions of pressure, density, velocity, total internal energy, Mach 
number, total pressure, and total temperature are presented in Figs. 28 
and 29. The key point to observe is the manner in which the terminal 
normal shock moves forward and amplifies in strength due to its move- 
ment. This is especially dramatic in the total pressure results of Fig. 
29, where the total pressure behind the moving normal shock reaches a 
level of 1.2 times the free-stream total pressure; in a similar manner 
the total temperature of Fig. 29 reaches a value of 1.2 times the 
free-stream total temperature. 

Another characteristic of inlet hammershock flow fields is the low 
subsonic Mach number (on the order of 0.2) behind the moving shock, 
as can be seen from Fig. 29; for comparison, the aftershock Mach 
number for a stationary normal shock at a Mach number of 2.50 is 
0.5130 and at an infinite Mach number is the limiting value 0.378. 

The hammershock calculation was performed using the implicit split 
characteristics solution algorithm with 41 grid points. Fig. 30 shows 
the Courant number distribution as a function of time. The inlet plane 
Courant number was reduced from a value slightly over 5 to 0.9 when 
the moving shock reached a x-location of approximately 3 in order to 
accurateiy track the moving shock upstream of the throat. 

5.4 40-60 INLET UNSTART/RESTART TRANSIENT 

During normal, started operation of mixed-compression inlets, 
variations in engine airflow are generally compensated for by the over- 
board bypass system. In this manner, the terminal shock can be 
positioned such that optimum inlet performance is maintained. If the 
bypass doors are scheduled to open upon unstart, the scheduled area 
must be large enough to compensate for the maximum reduction in 
engine airflow that could occur as a result of the unstart. However 
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FIGURE 28. FLOW VARIABLES DURIMG HAMMERSHOCK TRANSIENT 















FIGURE 3Q. COURANT NUMBER DISTRIBUTION DURING 
HANtMERSHOCK TRANSIENT 






scheduling bypass area for a bypass flow greater than that required 
would result in higher distortion and lower pressure recovery. If the 
engine continues to operate normally after the unstart transient, the 
increase in distortion could initiate a compressor stall. 

The present section will computationally examine the unstart/ 
restart transient operating characteristics of the 40-60 axisymmetric 
mixed-comprossion inlet with a cold pipe extension at a free-stream 
Mach number of 2.50. Experimental data reported in Ref. 29 will be 
used for comparison with calculated results. Details of the 40-60 inlet 
configuration and dynamic pressure transducer locations have been 
previously given in Fig. 22. 

Inlet unstart followed by controlled restart is shown in Fig. 31 
where the top two curves (disturbance door area and control door area) 
give the bypass door schedule in time, and the third curve from the 
top (centerbody position) represents centerbody translation in time from 
the fully retracted position. For purposes of the present one- 
dimensional computation, the disturbance door and control door area 
schedules were lumped together to provide a composite bypass door area 
schedule. The following is a description of the major events that 
occurred during the un start/ restart transient shown in Fig. 31: 

1. The terminal normal shock was initially placed on the verge of 
unstart, i.e., located immediately downstream of the inlet 
throat. 

2. A pulse decrease in disturbance door area caused the terminal 
normal shock to move upstream. 

3. The control doors are opened (increased area) to restore 
shock position in order to prevent inlet unstart. 

4. Despite the terminal shock control action, a net decrease in 
bypass mass flow was created and the inlet unstarted. 

5. The bypass doors were fully opened and the centerbody was 
translated forward to restart the inlet. 

6. Upon inlet restart, the centerbody was retracted to its ori- 
ginal position, and the control door area was decreased to 
move the terminal normal shock to a stable location down- 
stream of the inlet throat. 


DATA (Ref. 29) 
PROGRAM LAPIM 



FIGURE 31. IMLET UN START/RE START HISTORY 





In general, the present one-dimensional unsteady inviscid simulation 
agrees reasonably well with experiment. The disagreement shown in 
cowl lip pressure p^^ during unstart is due to viscous effects (flow 
separation caused by interaction of the unstarted normal shock with the 
centerbody boundary layer) influencing the "effective" area distribution 
of the inlet flow field upstream of the throat. The physical time for 
occurrence of both unstart and restart is accurately modeled by the 
simulation. After restart the inviscid simulation indicates some normal 
shock movement which is not supported by the experimental data. This 
movement is probably the result of viscous effects influencing shock 
location in the inlet diffuser through boundary-layer/shock interaction. 

Given that the present one-dimensional unsteady inviscid model is 
a reasonably accurate simulation of supersonic inlet un start/ restart phe- 
nomena as shown above, it is now in order to examine the computational 
details of the inlet flow field. Presented in Fig. 32 are the inlet 
area, static pressure, Mach number, and mass flow distributions at a 
time of 0.07 seconds which is just prior to initiation of bypass door 
closure. It should be noted that the downstream area distribution 
reflects struct blockage effects; however, bypass plenum effects are not 
included in the area distribution as all bypass flow phenomena are 
modeled as discussed earlier. Further note that the normal shock is 
positioned in the throat region and maintained there by throat bleed. 
The large effect of bypass mass removal on downstream mass flow is 
strikingly apparent in Fig. 32. 

Figure 33 shows the calculated details of the inlet unsfart tran- 
sient, which is extremely rapid (on the order of 0.01 to 0.u2 seconds). 
Note the decrease in inlet diffuser static pressure following unstart; 
this decrease in pressure is the direct result of the sizable decrease 
(50 to 70 percent) in inlet mass flow due to unstart mass spillage over 
the cowl. Further note the negative Mach number and mass flow within 
the inlet at a time of 0.0997 seconds. This is a direct result of the 
strong unstart transient creating a hammershock-like flow which has 
sufficient shock speed "drag" to reverse the local flow behind the 
moving shock. 
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Centerbody translation effects are shown in Fig. 34 where the 
forward movement serves to increase the throat area and decrease the 
cowl capture area. This physical geometry change results in the flow 
choking at the throat with supersonic diffuser flow terminating through 
a normal shock; note the rearward and then forward movement of this 
shock reflecting the changing area distribution due to centerbody 
translation. The mass flow essentially stabilizes throughout the inlet, 
although the mass flow entering the inlet decreases. 

Details of inlet restart are given Fig. 35. The actual restart is 
the result of decreasing the mass flow entering the inlet until such a 
value is reached that can be passed through the (choked) throat area. 
When this occurs the unstart normal shock moves rapidly into the inlet, 
passes through the throat, and creates a terminal normal shock in the 
inlet diffuser. 

With the inlet restarted, the centerbody is retracted as shown in 
Fig. 36 resulting in a decrease in throat area and an increase in 
capture area. This movement results in some fore and aft oscillations 
of the terminal normal shock as the mass flow through the inlet in- 
creases. At the last time shown (1.3951 seconds), the inlet is 
approaching a steady-state operating condition. 

Temporal details of the present un start/ restart calculation are 
shown in Figures 37 through 39 for four axial locations (3.22, 3.60, 
3.92, and 4.38) corresponding to locations where experimental pressure 
measurements were reported in Fig. 31. Note that these locations are 
in the throat region of the inlet under design conditions, i.e., no 
centerbody translation, as can be seen from reference to Fig. 32. With 
respect to the unstart/ restart transient, it is to be remembered from 
the previously discussed results that inlet unstart occurred at a time of 
approximately 0.09 seconds with inlet restart at a time of approximately 
0.65 seconds. 

Details of the inlet area variation with time are shown in Fig. 37. 
Forward/ rearward translation of the centerbody is apparent, with 
resulting increase/decrease of the throat area. 
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FIGURE 36.CENTERRODY RETRACTION BETWEEN TIME = a695t SECONDS 
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Fluid dynamic details of the transient un start/ re' cart flow are 
given in Figs. 38 and 39. The transient variation in all the fluid 
dynamic variables during unstart (at a time of approximately 
0.09 seconds) and during restart (at a time of approximately 
0.65 seconds) are quite pronounced. In addition, the transient oscil- 
lation of the terminal normal shock near the throat region after restart 
(time greater than approximately 0.65 seconds) can be seen quite 
clearly in Fig. 39 for the Mach number. 

The above calculation was performed using a prescribed corrected 
mass flow exit plane boundary condition. A Courant number condition 
of approximately ten (10) was utilized for all started phases of the cal 
culation with a Courant number of 0.90 utilized for all unstart and 
restart portions of the calculation where che normal shock moved 
rapidly. A total of 105 grid points was applied to discretize the flow 
field within the inlet inself; an additional seven (7) grid points were 
appended in front of the inlet for the un start/ restart phase. The 
"spike" which appears in the mass flow calculations presented earlier is 
due to the shock capturing nature of the numerical solution algorithm 
across strong normal shocks. A total of approximately 5,000 time steps 
was required to perform the complete transient calculation over a 1.40 
second physical time interval with a corresponding computer CPU time 
requirement of approximately 50 minutes on a Digital Equipment 
Corporation VAX-11/780. 

5.5 40-60 INLET DYNAMIC RESPONSE 

The experimental investigation of Ref. 32 was concerned with 
evaluation of the dynamic response characteristics of the 40-60 inlet 
configuration at a free-stream Mach number of 2.50. The inlet was 
coupled to a long cold pipe and subjected to sinusoidal bypass area 
disturbances as well as sinusoidal free-stream Mach number perturba- 
tions created via a gust plate. In the present section, frequency 
response parameters based upon Program LAPIN calculated results will 
be presented in comparison with the above discussed experimental 
results . 
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Figure 40 shows the comparison of calculated versus experimental 
results for the frequency response parameters of amplitude ratio and 
phase shift at various locations in the inlet and cold pipe. The cal- 
culations have been performed for three values of sinusoidal bypass 
area disturbance frequencies, namely 30, 50, and 100 hertz. Shock 
position amplitude ratio results indicate that the Program LAPIN shock 
movement is more heavily damped than experimentally observed; this is 
probably due to the neglect of viscous effects in the LAPIN inviscid- 
only treatment. A notable feature of the comparison is the extreme 
resonance in the 40-50 hertz range exhibited by both the experimental 
data and the inviscid calculations. 

Calculated frequency response characteristics of the 40-60 inlet, 
with long cold pipe to upstream Mach number sinusoidal oscillations are 
presented in Fig. 41. The calculations have been performed at a 
frequency of 10 hertz. Calculated pressure response at various loca- 
tions in the inlet are in excellent agreement with experiment, both for 
magnitude and phase. The experimental shock damping (amplitude 
ratio) is once again underpredicted. 

5.6 60-40 INLET 

The NASA LeRC 60-40 inlet is an axisymmetric mixed-compression 
inlet with 60 percent of the effective supersonic area contraction 
occurring externally and 40 percent supersonic area contraction occurr- 
ing internally at the design free-stream Mach number of 2.50. The 
external supersonic compression is accomplished via a biconic center- 
body configuration. Details of the 60-40 inlet configuration are given 
in Fig. 42. The inlet was equipped with a translating centerbody, 
throat region cowl and centerbody bleed, and overboard bypass doors. 
Pertinent details concerning experimental investigations of this inlet 
conducted in the NASA LeRC 10- by 10-foot Supersonic Wind Tunnel 
are contained in Ref. 34 through 37. 
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5.7 60- 0 INLET STEADY-STATE PERFORMANCE 


The present sectiofi will computationally examine the steady-state 
operation of the 60-40 axisymmetric mixed-compression inlet at a free- 
stream Mach number of 2.50. Experimental data for a porous bleed 
configuration reported in Ref. 37 will be used for comparison with 
calculated results. Details of the 60-40 inlet configuration have been 
previously given in Fig. 42. 

Overall inlet performance in terms of total pressure recovery as a 
function of engine mass flow variation is shown in Fig. 43. The calcu- 
lated total pressure recovery is in good agreement with the experimental 
data, both in distribution and magnitude. Also shown in Fig. 43 is the 
bleed plenum total pressure recovery as a function of throat bypass 
bleed mass flow ratio. The linear characteristics of bleed plenum total 
pressure recovery with respect to variation in throat bypass bleed mass 
flow is well predicted by the present inviscid bleed/bypass model. 

Inlet static pressure distributions for various mass flow ratios are 
given in Fig. 44 for both the internal cowl and centerbody surfaces. 
As with the 40-60 inlet discussed previously, the calculated terminal 
normal shock is significantly downstream of the experimentally indicated 
location. As noted earlier, this discrepancy can be traced to the 
inviscid treatment of mass bleed effects without any allowance for 
viscous (boundary-layer) interaction. Irrespective of this, the calcu- 
lated static pressure distribution is in excellent agreement with 
experiment in regions upstream of the bleed locations and downstream of 
the calculated terminal shock. 

5.8 60-40 INLET TRANSIENT UNSTART LIMITS 

Program LAPIN calculations of transient unstart limits for the 60-40 
inlet with porous throat bypass were computed and compared to experi- 
mental data from Ref. 37. These results are shown in Fig. 45 for two 
bleed plenum volumes. It should be noted that the transient pulse was 
generated differently between experiment and simulation. During the 
experiment the bypass door opening was pulsed and then related to 
Inlet corrected airflow through a steady-state correlation to obtain the 
stability index. In the case of the simulation, the corrected airflow 
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W£)s pulsed directly, Steady“state and low plenum volume results are in 
good agreement with experiment. The large plenum simulation is some- 
what lower than experiment, with the difference suggesting that the 
transient pulse generation technique becomes more of a factor at the 
high airflow variations. 

5.9 B-70 INLET 

The B-70 inlet is a two-dimensional, external-internal compression 
inlet with bypass doors, variable throat area, throat bleed, and three 
external ramps. Details of the B-70 inlet are given in Fig. 46. Per- 
tinent details concerning experimental investigation of this inlet in the 
AEDC 16-Foot Supersonic Wind Tunnel are given in Ref. 38. 

5.10 B-70 INLET STEADY-STATE PERFORMANCE 

The present section will computationally examine the steady-state 
operation of the B-70 two-dimensional inlet at a free-stream Mach 
number of 3.0. Experimental data reported in Ref. 38 will be used for 
comparison with calculated results. Details of the B-70 inlet con- 
figuration have been previously given in Fig. 46. 

Figure 47 presents a comparison of the calculated steady-state 
static pressure distribution relative to experimental measurements. As 
can be seen, the experimental terminal normal shock is located in the 
throat region while the calculated normal shock is located downstream of 
the throat. In addition, the experimental pressure level is some 7% 
higher than the calculation downstream of the normal shock. These 
results are related in trend with the axisymmetric 40-60 inlet results 
reported earlier in this report and reflect the strong influence of 
viscous (boundary-layer) interaction effects in the throat region and 
downstream which are not accounted for in the present inviscid 
calculation . 
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FIGURE 46. B-70 INLET GEOMETRY 













6.0 CONCLUDING REMARKS 


The results presented in this report indicate that the transient 
flow field in a supersonic mixed-compression inlet with large flow per- 
turbations (hammershock, un start/ restart, etc.) can be reasonably well 
simulated via a quasi-one-dimensional inviscid approach using a shock 
capturing numerical solution algorithm. In order to further improve the 
applicability and accuracy of this approach, viscous boundary-layer 
effects and their interaction with the inviscid flow must be properly 
incorporated into the analysis. 

Five different numerical algorithms have been utilized in the 
present work with a summary given in Table 2. For supersonic inlet 
applications with embedded moving shock waves, the split characteris- 
tics approach provides the highest fidelity numerical solution and is 
recommended for practical applications. 
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APPENDIX A 

Straight Sonic Line Correction Factor Derivation 


Consider a two-dimensional (j=0) or axisymmetric (j=1) flat-faced 
body of half-height/ radius a in a supersonic flow as shown below. 


•SHOCK 


. ASSUMED STRAIGHT 
SONIC LINE (* CONDITIONS) 


Under the assumption that the sonic line is straight as illustrated 
above, continuity considerations require that 


S* __ 1 P OO V< 
"i" " 2 ! p*V* 


(A-1) 


where the * indicates evaluation at sonic (M=1.0) conditions. The 
rations p^/p* and V^/V* are easily determined from classical 
gasdynamics . 

However, the actual sonic line is not straight but curved as shown 
below ^^SHOCK 

^CURVED 

/ X SONIC UNE 


L-.|. 


4 





Vv. A 1 


OtAs.,..- .. . . 

OF H-'A ' ■ 

Thus the actual flow crossing the dashed line is, in reality, not at 
sonic conditions and thus 

i_ - 2 - (A-2) 

• 2» pV 


where the mass flux pV crossing the dashed line is evaluated at the 
corresponding local conditions. Forming the ratio of Eq. (A-1) 
and (A-2) yields 

pV »Fp*V (A-3) 


where 



can be considered as a straight sonic line correction factor which, when 
multipled by the corresponding straight sonic line mass flux p*V*, 
yields the actual mass flux pV crossing the dashed line shown above. 

For a two-dimensional or axisymmetric flat-faced body in super- 
sonic flow, the actual shock wave standoff distance 5 can be easily 
determined following Moeckel (Ref. 3) while the straight sonic line 
standoff distance 6* is given by Eq. (A-1). In this manner the 
straight sonic line correction factor F given by Eq. (A-4) becomes a 
function of free-stream Mach number only and is given in Fig. A-1 for 
both two-dimensional and axisymmetric bodies. 



PIOURE A-1. STRAIGHT SONIC LINE CORRECTION FACTOR 
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For the case of mass flow p|U| Into the flat-faced body as shown 
below 



continuity considerations for the assumed straight sonic line require that 


S* 1 P oo V oo 

a ” 2J p*V* 


P|U| 

P oo V oo 


(A-5) 


In a similar manner for the actual flow crossing the dashed line 


^ 1 PooVoo 

a 2i pV 


P|U| 

P oo V oo 


(A-6) 


Forming the ratio of Eq. (A-5) and (A-6) yields exactly the same 
results as Eq. (A-3) and (A-4) presented previously. The 
important point here is that the straight sonic line correction factor 
F is unaffected by mass flow through the body. 
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APPENDIX B 
SYMBOLS 

Cross-sectional area; tridiagonal matrix coefficient 

Speed of sound; body radius, geometry parameter in Fig. 7 

Tridiagonal matrix coefficient 

Tridiagonal matrix coefficient 

Specific heat at constant pressure 

Specific heat at constant volume 

Diagonal matrix 

Diagonal components of diagonal matrix 
Total interna! energy 
Specific internal energy 
Flux vector 

Momentum equation source term 

Right-hand-side vector 

Total enthalpy 

Identity matrix 

Jacobian matrix 

Flow coefficient 

Length 

Mach number 

Continuity equation source term 
Pressure 

Characteristic component; equivalent one-dimensional flux 
quantities 

Energy equation source term 
Gas constant 
Right-hand-side vector 
Diagonalizing matrix 
Temperature 
Time 

Solution vector 
Primitive solution vector 
Reference velocity 
Axial velocity 
Velocity; normal velocity 
volume 
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APPENDIX B 
SYMBOLS 
(continued) 

Axial distance aiong centeriine 
Normai distance from centerline 
Cowl radius 
Specific heat ratio 
Shock stand-off distance 
Switching parameter 
Local flow angle 
Diagonalized Jacobian matrix 
Eigenvalue 

Nondimensional axial distance 
Density 

Nondimensional time 

Decomposition function; conical flow angle 


Subscript 


c 

Cp 

L 

P 

ref 

s 

T, O 
1,2 
1,2,3 


Cone surface 

Conservative to primitive variables 

Local 

Plenum 

Reference 

Shock 

Total 

State ahead of and behind shock (scalar) 
Mass, momentum, energy (vector) 
Free-stream 


Superscript 

* Sonic; known co.istraint value 

* Positive characteristic value 
Negative characteristic value 

" Dimensional variable; Runge-Kutta predictor value 

A Nondimensional variable; factored solution variable 

Vector quantity 
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APPENDIX 3 
SYMBOLS 
(concluded) 
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Mathematical Operators 
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6“ 


Backward difference spatiai operator 
Forward difference spatiai operator 
Fourth-order dissipation difference operator 
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